library(tidyverse)
library(tidymodels)
library(modelr)
library(ISLR)
library(GGally)
library(pROC)
library(cowplot)
library(OneR)
library(rlang)
library(caret)
set.seed(1992)
La regresión logística es útil para problemas de predicción de clases. En este caso vamos a tratar de resolver el problema de predecir si una persona va defaultear su deuda de tarjeta de crédito en base a ciertos predictores.
Este conjunto de datos proviene de la librería ISLR (Introduction to Statistical Learning Using R) de James, Witten, Hastie y Tibshirani.
# cargamos los datos y observamos su estructura
default <- Default
glimpse(default)
Rows: 10,000
Columns: 4
$ default <fct> No, No, No, No, No, No, No, No, No, No, No, No, No, N...
$ student <fct> No, Yes, No, No, No, Yes, No, Yes, No, No, Yes, Yes, ...
$ balance <dbl> 729.5265, 817.1804, 1073.5492, 529.2506, 785.6559, 91...
$ income <dbl> 44361.625, 12106.135, 31767.139, 35704.494, 38463.496...
Tiene 4 variables:
default: La clase que queremos predecir.
student: Binaria que indica si la persona es estudiante.
balance: Balance promedio que le queda a la persona luego de sus pagos mensuales.
income: Ingreso de la persona.
Analicemos la distribución de la clase.
default %>%
group_by(default) %>%
summarise(numero_casos=n())
`summarise()` ungrouping output (override with `.groups` argument)
Vemos que estamos trabajando con un problema de clasificación con un claro desbalance de clase.
Realizamos un gráfico exploratorio completo para ver el comportamiento y las relaciones entre las variables. El color rojo designa a quienes no defaultean y el azul a los que sí.
# graficamos con ggpairs coloreando por variable a predecir
g <- ggpairs(default, 1:4, title = "Correlograma de variables",
mapping = aes(colour= default)) +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
theme_bw()
# hacemos un loop para cambiar los colores del gráfico
for(i in 1:g$nrow) {
for(j in 1:g$ncol){
g[i,j] <- g[i,j] +
scale_fill_brewer(palette="Set1") +
scale_color_brewer(palette="Set1")
}
}
g
¿Qué pueden decir de la relación entre balance y default?
¿Y entre income y balance?
¿Cuáles parecen ser buenas variables para discriminar entre quienes defaultean y quienes no?
Para modelizar va a ser necesario tener la variable default como numérica. Definimos la variable como {0,1} para los valores {“No”,“Yes”}
default <- default %>%
mutate(default = case_when(default == "No" ~ 0,
default == "Yes" ~ 1))
default
Queremos estimar \(P(Default=Yes|X)=P(X)\) para cada individuo y a partir de ello poder definir un punto de corte para predecir quiénes son los que van a entrar en default.
En este caso estamos modelando la probabilidad de la siguiente manera:
\(P(X)= \beta_0 + \sum\limits_{j=1}^p \beta_j X\)
Veamos que tan bueno es el modelo lineal para esto, usando balance como predictor.
test_mco <- default %>%
lm(formula = default ~ balance, data = .)
tdy = test_mco %>% tidy()
tdy
test_mco %>% glance()
Ambos estimadores son significativos y el test de significatividad global del modelo también es significativo.
Veamos un gráfico de nuestro modelo.
Parece tener bastantes problemas para estimar la probabilidad de default de los individuos. Por ejemplo, vemos que hay varios individuos a los cuales les asigna una probabilidad negativa.
Para evitar estos problemas, usamos la función logística.
\(P(Y=1|X)= \frac{e^{\beta_0 + \sum\limits_{j=1}^p \beta_j X}}{1+e^{\beta_0 + \sum\limits_{j=1}^p \beta_j X}}\)
El lado derecho se llama expit
Esta función acota el resultado entre 0 y 1, lo cual es mucho más adecuado para modelar una probabilidad.
Luego de hacer algunas operaciones, podemos llegar a la expresión:
\(\log {\frac{P(x)}{1-P(x)}}= \beta_0 + \sum\limits_{j=1}^p \beta_j X\)
El lado izquierdo es el logaritmo de los odds y se llama logit.
Realizamos una partición una partición entre dataset de entrenamiento (70%) y testeo (30%) usando la función initial_split del paquete rsample de tidymodels.
# fijamos semilla
set.seed(44)
# Partición Train y Test, indicando proporción
train_test <- initial_split(default, prop = 0.7)
# armamos dataframe de testeo y entrenamiento
default <- training(train_test)
test <- testing(train_test)
# vemos el contenido
default %>%
dim_desc() # 7000 filas
[1] "[7,000 x 4]"
test %>%
dim_desc() # 3000 filas
[1] "[3,000 x 4]"
Para aplicar la regresión logística primero usamos la función formulas del paquete modelr para crear un objeto que contiene todas las fórmulas que vamos a utilizar. En .response especificamos la variable respuesta de nuestras fórmulas y luego nombramos las fórmulas que queramos armar.
Así, armaremos distintos modelos combinando distintas variables.
logit_formulas <- formulas(.response = ~ default, # único lado derecho de las formulas.
bal = ~ balance,
stud = ~ student,
inc = ~ income,
bal_stud = ~ balance + student,
bal_inc = ~ balance + income,
stud_inc = ~ student + income,
full = ~ balance + income + student
)
logit_formulas # observamos el objeto formulas
$bal
default ~ balance
$stud
default ~ student
$inc
default ~ income
$bal_stud
default ~ balance + student
$bal_inc
default ~ balance + income
$stud_inc
default ~ student + income
$full
default ~ balance + income + student
Procedemos a crear los modelos a partir de estas fórmulas.
models <- data_frame(logit_formulas) %>% # dataframe a partir del objeto formulas
mutate(models = names(logit_formulas), # columna con los nombres de las formulas
expression = paste(logit_formulas), # columna con las expresiones de las formulas
mod = map(logit_formulas, ~glm(., family = 'binomial', data = default))) # Que estamos haciendo acá? Que vamos a encontrar en la columna?
models
La funcíón glm() nos permite crear un modelo lineal generalizado (Generalized Linear Model). Al igual que la función lm() toma como argumentos una formula y los datos pero también se debe especificar el argumento family: indicamos la distribución del error y la función link que vamos a utilizar en el modelo.
Algunas familias son:
Binomial: link=logit
Poisson: link=log
Gaussiana: link=identidad
Como estamos trabajando con un fenómeno que suponemos tiene una distribución binomial, así lo especificamos en el parámetro family.
Probamos los primeros tres modelos, aquellos que tienen un único predictor. Usamos la función tidy para obtener los parámetros estimados para estos tres modelos.
models %>%
filter(models %in% c('bal','stud','inc')) %>%
mutate(tidy = map(mod, tidy)) # Qué realizamos en este paso? Que va a tener esta columna?
Para acceder a los elementos de la nueva columna tidy debemos desanidarla.
models %>%
filter(models %in% c('bal','stud','inc')) %>%
mutate(tidy = map(mod, tidy)) %>% # Qué realizamos en este paso? Que va a tener esta columna?
unnest(tidy) %>%
mutate(estimate=round(estimate,5), # redondeamos valores para facilitar lectura
p.value=round(p.value,4))
Observamos que todos los modelos tienen coeficientes significativos excepto el coeficiente asociado al Ingreso (income), el cual tiene un p-valor de 0.17.
Recordando la ecuación para modelar la probabilidad:
\(P(Y=1|X)= \frac{e^{\beta_0 + \sum\limits_{j=1}^p \beta_j X}}{1+e^{\beta_0 + \sum\limits_{j=1}^p \beta_j X}}\)
Se observa que ahora las variables ya no tienen una relación lineal con la probabilidad. En este modelo un coeficiente positivo indica que frente a aumentos de dicha variable la probabilidad aumenta, mientras que un coeficiente negativo nos indica lo contrario. Para nuestro caso:
El coeficiente de 0.0056 de la variable balance indica que cuando el balance aumenta la probabilidad de default aumenta.
El coeficiente de 0.40886 de la variable studentYes indica que cuando una persona es estudiante su probabilidad de default aumenta respecto a una persona que no lo es.
El coeficiente de -0.00001 de la variable income indica que cuando el ingreso aumenta la probabilidad de default disminuye.
Recomendamos leer el capítulo de regresión logística de Interpretable Machine Learning: A Guide for Making Black Box Models Explainable de Molnar Cristoph para una discusión más profunda de la interpretación de un modelo de regresión logística.
Ahora probamos con un modelo que utiliza las tres predictoras.
models %>%
filter(models == "full") %>%
mutate(tidy = map(mod,tidy)) %>%
unnest(tidy) %>%
mutate(estimate=round(estimate,5),
p.value=round(p.value,4))
En el modelo con las tres variables los coeficientes de studentYes e income han cambiado de signo, por lo cual su efecto sobre la probabilidad de default es inverso al que observamos en los modelos simples. También hay que agregar ahora que el efecto de la variable se observa al estar controlando por las dos covariables restantes.
Asimismo, se observa que el coeficiente de income no resulta significativo tampoco en este modelo (p-valor>0.05).
Con map() agregamos la función glance para traer información relevante para la evaluación del modelo. Con unnest() accedemos a dicha información. Por último, agregamos una columna con el porcentaje de deviance explicado por cada modelo y ordenamos el dataset según su valor de deviance.
# Calcular las medidas de evaluación para cada modelo
models <- models %>%
mutate(glance = map(mod,glance))
# Obtener las medidas de evaluacion de interes
models %>%
unnest(glance) %>%
# Calculamos la deviance explicada
mutate(perc_explained_dev = 1-deviance/null.deviance) %>%
select(-c(models, df.null, AIC, BIC)) %>%
arrange(deviance)
El modelo que utiliza las 3 variables es el que minimiza el deviance, aunque el modelo que excluye a la variable income (poco significativa) está muy cerca. Los 3 últimos modelos reducen muy poco el deviance respecto a la deviance nula.
Realizamos los gráficos para el modelo completo y uno de los modelos con mayor deviance (student+income).
Comenzamos agregando las predicciones con augment con el parámetro type="response". La función augment hereda el argumento type.predict de la función predict.
Si type.predict = 'link' la predicción es en términos de la función link. En nuestro caso son el logaritmo de las odds, es decir, los valores que toma la expresión logit.
Si type.predict = 'response' la predicción son las probabilidades de que la observación pertenezca a la clase positiva. En nuestro caso, devuelve la probabilidad de la que persona defaultee.
# Añadir las predicciones
models <- models %>%
mutate(pred= map(mod,augment, type.predict = "response"))
#Observaciones con probabilidad más baja
models$pred[1]$bal %>% arrange(.fitted) %>% head(10)
#Observaciones con probabilidad más alta
models$pred[1]$bal %>% arrange(desc(.fitted)) %>% head(10)
Guardamos las predicciones para los modelos mencionados.
# Modelo completo
prediction_full <- models %>%
filter(models=="full") %>%
unnest(pred, .drop=TRUE)
#Modelo malo
prediction_bad <- models %>%
filter(models=="stud_inc") %>%
unnest(pred, .drop=TRUE)
# graficamos el modelo completo
violin_full = ggplot(prediction_full, aes(x=default, y=.fitted, group=default, fill=factor(default))) +
geom_violin() +
theme_bw() +
guides(fill=FALSE) +
labs(title='Violin plot', subtitle='Modelo completo', y='Predicted probability')
# graficamos el modelo malo
violin_bad = ggplot(prediction_bad, aes(x=default, y=.fitted, group=default, fill=factor(default))) +
geom_violin() +
theme_bw() +
guides(fill=FALSE) +
labs(title='Violin plot', subtitle='Modelo malo', y='Predicted probability')
# mostramos ambos
plot_grid(violin_bad, violin_full)
En los gráficos de violin observamos:
En el eje de abscisas la clase verdadera: Default o No Default (1 o 0).
En el eje de ordenadas la probabilidad predicha por nuestro modelo.
El gráfico nos muestra la distribución de la cantidad de observaciones por su clase real y la probabilidad que le asigna nuestro modelo.
¿Cuál parece ser un punto de corte adecuado para cada modelo?
Se genera una función para realizar un gráfico de Hosmer-Lemeshow para un dataset. Para ello se fijan los siguientes parámetros:
dataset: conjunto de datos
predicted_column: columna con la probabilidad predicha
class_column: columna con la clase a predecir
possitive_value: valor de la clase a predecir
bins: cantidad de grupos del gráfico
color: color de los puntos
nudge_x: desplazamiento de la etiqueta en el eje x
nudge_y: desplazamiento de la etiqueta en el eje y
Hosmer_Lemeshow_plot <- function(dataset, predicted_column, class_column, bins, positive_value, color='forestgreen', nudge_x=0, nudge_y=0.05){
# Asignar los grupos a las observaciones de acuerdo a la probabilidad predicha
dataset['group'] <- bin(dataset[predicted_column], nbins = bins, method = 'l', labels=c(1:bins))
# Contar la cantidad de casos positivos por grupo
positive_class <- dataset %>% filter(!!sym(class_column)==positive_value) %>% group_by(group) %>% count()
# Obtener la media de las predicciones por grupo
HL_df <- dataset %>% group_by(group) %>% summarise(pred=mean(!!sym(predicted_column)), count=n()) %>%
inner_join(.,positive_class) %>%
mutate(freq=n/count)
# Gráfico
HM_plot <- ggplot(HL_df, aes(x=pred, y=freq)) +
geom_point(aes(size=n), color=color) +
geom_text(aes(label=n),nudge_y = nudge_y)+
geom_abline(slope = 1, intercept = 0, linetype='dashed') +
theme_bw() +
labs(title='Hosmer-Lemeshow', size='Casos', x="Probabilidad Predicha", y="Frecuencia observada")
return(HM_plot)
}
Generamos los gráficos pasandole lo parámetros.
# modelo completo
Hosmer_Lemeshow_plot(prediction_full, '.fitted', 'default', 10, 1) +
labs(subtitle="Modelo completo")
# modelo malo
Hosmer_Lemeshow_plot(prediction_bad, '.fitted', 'default', 10, 1, color = "firebrick", nudge_y = 0.003) + scale_x_continuous(limits = c(0.02,.06)) + scale_y_continuous(limits = c(.02,.06)) + labs(subtitle="Modelo malo")
En los gráficos de Hosmer-Lemeshow observamos:
En el eje de abscisas la probabilidad predicha de default.
En el eje de ordenadas la frecuencia de clase, el cociente entre cantidad de individuos que están en default y el total de individuos.
La línea punteada designa la igualdad entre probabilidad predicha y frecuencia de clase.
Los círculos, que se construyen de la siguiente manera:
Aquellos círculos que se ubiquen por encima de la línea punteada indican que el modelo está subestimando la probabilidad para dichos grupos. Mientras que si los círculos se ubican por debajo el modelo está sobreestimando la probabilidad para dichos grupos.
¿Para qué valores parece existir una sobreestimación de la probabilidad? ¿Para cuáles subestimación?
# Calculamos curvas ROC
roc_full <- roc(response=prediction_full$default, predictor=prediction_full$.fitted)
Setting levels: control = 0, case = 1
Setting direction: controls < cases
roc_bad <- roc(response=prediction_bad$default, predictor=prediction_bad$.fitted)
Setting levels: control = 0, case = 1
Setting direction: controls < cases
Graficamos ambas en un mismo plot.
ggroc(list(full=roc_full, bad=roc_bad), size=1) +
geom_abline(slope = 1, intercept = 1, linetype='dashed') +
theme_bw() +
labs(title='Curvas ROC', color='Modelo')
print(paste('AUC: Modelo completo', round(roc_full$auc,3)))
[1] "AUC: Modelo completo 0.952"
print(paste('AUC: Modelo malo', round(roc_bad$auc,3)))
[1] "AUC: Modelo malo 0.562"
¿Qué significa cada uno de los ejes?
Hasta ahora hemos evaluado el modelo de manera general, pero el resultado final del modelo debe consistir en asignar a la persona una clase predicha. En nuestro caso debemos establecer un punto de corte según el cual vamos a separar a las personas en quienes defaultean y quienes no.
Probamos varios puntos de corte y graficamos el accuracy, la sensibilidad o recall, la especificidad y la precisión para cada uno de ellos.
| Clases predichas / Clases | Positiva | Negativa |
|---|---|---|
| Positiva | True Pos | False Pos |
| Negativa | False Neg | True Neg |
Recordemos que:
\(accuracy = \frac{TP+TN}{TP+FP+FN+TN}\)
\(sensitivity = recall = \frac{TP}{TP+FN}\)
\(specificity = \frac{TN}{TN+FP}\)
\(precision = \frac{TP}{TP+FP}\)
prediction_metrics <- function(cutoff, predictions=prediction_full){
table <- predictions %>%
mutate(predicted_class=if_else(.fitted>cutoff, 1, 0) %>% as.factor(),
default= factor(default))
confusionMatrix(table(table$predicted_class, table$default), positive = "1") %>%
tidy() %>%
select(term, estimate) %>%
filter(term %in% c('accuracy', 'sensitivity', 'specificity', 'precision')) %>%
mutate(cutoff=cutoff)
}
cutoffs = seq(0.01,0.95,0.01)
logit_pred= map_dfr(cutoffs, prediction_metrics)%>% mutate(term=as.factor(term))
ggplot(logit_pred, aes(cutoff,estimate, group=term, color=term)) + geom_line(size=1) +
theme_bw() +
labs(title= 'Accuracy, Sensitivity, Specificity y Precision', subtitle= 'Modelo completo', color="")
¿Qué podemos observar en el gráfico?
¿Podemos definir un buen punto de corte? ¿Cuál sería?
¿Por qué la especificidad tiene ese comportamiento?
Seleccionamos el modelo completo, ya que es el que maximizaba el porcentaje de deviance explicada y en base a lo que vimos definimos un punto de corte en 0.25 (pueden probar otros)
sel_cutoff = 0.27
# Creamos el modelo
full_model <- glm(logit_formulas$full, family = 'binomial', data = default)
# Agregamos la predicciones al dataset de testeo
table= augment(x=full_model, newdata=test, type.predict='response')
# Clasificamos utilizamos el punto de corte
table=table %>%
mutate(predicted_class=if_else(.fitted>sel_cutoff, 1, 0) %>% as.factor(),
default= factor(default))
# Creamos la matriz de confusión
confusionMatrix(table(table$default, table$predicted_class), positive = "1")
Confusion Matrix and Statistics
0 1
0 2840 55
1 51 54
Accuracy : 0.9647
95% CI : (0.9574, 0.971)
No Information Rate : 0.9637
P-Value [Acc > NIR] : 0.4092
Kappa : 0.4864
Mcnemar's Test P-Value : 0.7708
Sensitivity : 0.49541
Specificity : 0.98236
Pos Pred Value : 0.51429
Neg Pred Value : 0.98100
Prevalence : 0.03633
Detection Rate : 0.01800
Detection Prevalence : 0.03500
Balanced Accuracy : 0.73889
'Positive' Class : 1
Al explorar el dataset vimos que existía un fuerte desbalance de clase. Sólo el 3% de las observaciones pertenecen a personas que defaultearon. Esto puede tener un efecto en las estimaciones del modelo y su clasificación final.
Existen dos maneras sencillas con las cuales podemos trabajar con una clase desbalanceada:
Sobre-muestreo (oversampling) de la clase minoritaria
Sub-muestreo (undersampling) de la clase mayoritaria
La función glm puede tomar como argumento una columna (weigths) de ponderadores para poder hacer esto. Podemos asignar pesos mayores a 1 a la clase minoritaria (oversampling) o menores a 1 a la clase mayoritaria (undersampling). En nuestro problema vamos a realizar un sobresampleo de la clase minoritaria.
# Creamos la columna de ponderadores
default <- default %>% mutate(wt= if_else(default==1,20,1))
# Creamos los modelos con la data 'balanceada'
balanced_models <- data_frame(logit_formulas) %>% # dataframe a partir del objeto formulas
mutate(models = names(logit_formulas), # columna con los nombres de las formulas
expression = paste(logit_formulas), # columna con las expresiones de las formulas
mod = map(logit_formulas, ~glm(.,family = 'binomial', data = default, weights = wt))) #Pasamos la columna wt como ponderadores
Vemos las estimaciones de los parámetros para el modelo completo. ¿Existen cambios?
Ahora veamos la evaluación de los modelos ¿Qué pasó con el porcentaje de deviance explicada? ¿Y con la nula?